/*
 * slipmodel.c
 *
 *  Created on: Jul 23, 2013
 *      Author: bach
 *
 *  -------------------------------------------
 *  int readSlipmodel(struct slipmodel *pscmp);
 *  -------------------------------------------
 *  Reading a pscmp input file to a slipmodel struct, using filename stored in (struct) slipmodel.filename
 *  Automatically setting patch vector indices.
 *
 *
 *
 */


#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include "slipmodel.h"
#include "nrutil.h"
#include "c_recipes.h"
#include "dft2d.h"
#include "mathfunctions.h"

/******************************************************************/
/*                READ SLIPMODEL FROM FILE                        */
/******************************************************************/

int slipmodelRead(struct slipmodel *pscmp) {

    FILE *fin;
    char cdum[256], *cerror, *cdum2, *cdum3;
    int obs_type, no_outputs, segment, patches = 0, sum_no_patches = 0, i, j;
    fpos_t pos;


    // open input file
    fin = fopen((*pscmp).filename, "r");
    if (fin == NULL) {
        fprintf(stderr, "\nError opening slipmodel file %s\n(in function readSlipmodel)\n", (*pscmp).filename);
        return(1);
    }



    // Kommentarzeilen ueberspringen
    cerror = fgets(cdum, 256, fin);
    while (cdum[0] == '#')
        cerror = fgets(cdum, 256, fin);

    if (cerror == NULL)
            fprintf(stderr, "Error reading from pscmpfile!\n");

    sscanf(cdum, "%d", &obs_type); // observation type
    if (obs_type != 2) {
        printf("\nError: pscmpfile not using regular observation array!\nChange this source code to fit.\n");
        return(2);
    }

    // read two rows of the observation array
    cerror = fgets(cdum, 256, fin);
    sscanf(cdum, "%d %lf %lf", &(*pscmp).numberObsLat, &(*pscmp).obsLatMin, &(*pscmp).obsLatMax);
    cerror = fgets(cdum, 256, fin);
    sscanf(cdum, "%d %lf %lf", &(*pscmp).numberObsLon, &(*pscmp).obsLonMin, &(*pscmp).obsLonMax);

    // next row
    cerror = fgets(cdum, 256, fin);
    while (cdum[0] == '#')
        cerror = fgets(cdum, 256, fin);

    // prestress
    cerror = fgets(cdum, 256, fin);
    sscanf(cdum, "%*d %lf %lf %lf %lf %lf %lG %lG %lG", &(*pscmp).friction, &(*pscmp).skempton, &(*pscmp).strike, &(*pscmp).dip, &(*pscmp).rake, &(*pscmp).sigma1, &(*pscmp).sigma2, &(*pscmp).sigma3);

    cerror = fgets(cdum, 256, fin);

    cdum2 = strchr(cdum, '\'')+1;
    cdum3 = strrchr(cdum, '\'');
    cdum3[0] = '\0';
    strcpy((*pscmp).cmbOutFolder, cdum2);



    // OUTPUT Zeilen ueberspringen
    for (i = 1; i <= 7; i++)
        cerror = fgets(cdum, 256, fin);
    sscanf(cdum, "%d", &no_outputs);

    // Anzahl ausgelesener Outputfiles Zeilen ueberspringen
    //for (i = 1; i <= no_outputs; i++)
    cerror = fgets(cdum, 256, fin);
    cdum2 = strchr(cdum, '\'')+1;   // find first ' sign
    cdum3 = strrchr(cdum, '\'');    // find second ' sign
    cdum3[0] = '\0';                // replace second ' by string-end sign
    strcpy((*pscmp).cmbOutFile, cdum2);         // copy file name until newly added string-end sign


    cerror = fgets(cdum, 256, fin);

    while (cdum[0] == '#')
        cerror = fgets(cdum, 256, fin);

    // OUTPUTS (4 rows)
    for (i = 1; i <= 4; i++)
        cerror = fgets(cdum, 256, fin);

    while (cdum[0] == '#')
        cerror = fgets(cdum, 256, fin);

    // number of slip segments
    sscanf(cdum, "%d", &(*pscmp).numberSegments);

    cerror = fgets(cdum, 256, fin);
    while (cdum[0] == '#') {
        fgetpos(fin, &pos);
        cerror = fgets(cdum, 256, fin);
    }
    // Position zurueck setzen auf Zeilenanfang
    fsetpos(fin, &pos);

    // Vektor fuer Segmentdaten, pro Zeile ein Segment
    (*pscmp).segmentData = dmatrix(1, (*pscmp).numberSegments, 1, 10);
    (*pscmp).numberPatches = ivector(1, (*pscmp).numberSegments);



    /************************************************************************/
    /*                                                                      */
    /*                  Read Segment and Patch data                         */
    /*                                                                      */
    /************************************************************************/

    // read Segment data and determine amount of patches
    for (i = 1; i <= (*pscmp).numberSegments; i++) {
        // Segment Spezifikationen
        // lat   lon    depth length width strike dip   np_st np_di start_time
        // [deg] [deg]  [km]  [km]   [km]  [deg]  [deg] [-]   [-]   [day]
        cerror = fgets(cdum, 256, fin);
        sscanf(cdum, "%d %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf",
                &segment, &((*pscmp).segmentData)[i][1], &((*pscmp).segmentData)[i][2], &((*pscmp).segmentData)[i][3],
                &((*pscmp).segmentData)[i][4], &((*pscmp).segmentData)[i][5], &((*pscmp).segmentData)[i][6],
                &((*pscmp).segmentData)[i][7], &((*pscmp).segmentData)[i][8], &((*pscmp).segmentData)[i][9], &((*pscmp).segmentData)[i][10]);
        ((*pscmp).numberPatches)[i] = (int) (((*pscmp).segmentData)[i][8] * ((*pscmp).segmentData)[i][9]);
        sum_no_patches += ((*pscmp).numberPatches)[i];

        for (j = 1; j <= ((*pscmp).numberPatches)[i]; j++) {
            cerror = fgets(cdum, 256, fin);
        }
    }

    // Matrix fuer Patch data erzeugen
    (*pscmp).patchData = dmatrix(1, sum_no_patches, 1, 5);
    // File Pointer zuruecksetzen auf Beginn der Daten
    fsetpos(fin, &pos);
    patches = 0;

    for (i = 1; i <= (*pscmp).numberSegments; i++) {
        // Segment Daten ueberspringen, da bereits eingelesen
        cerror = fgets(cdum, 256, fin);

        // Patch Spezifikationen
        // pos_s   pos_d    slp_stk slp_dip open
        // [km]    [km]     [m]     [m]     [m]
        for (j = 1; j <= ((*pscmp).numberPatches)[i]; j++) {
            patches++;
            cerror = fgets(cdum, 256, fin);
            sscanf(cdum, "%lf %lf %lf %lf %lf", &((*pscmp).patchData)[patches][1], &((*pscmp).patchData)[patches][2], &((*pscmp).patchData)[patches][3], &((*pscmp).patchData)[patches][4], &((*pscmp).patchData)[patches][5]);
        }
    }





    if (sum_no_patches != patches) {
        printf("\nAnzahl Patches nicht konsistent!\n\n");
        return(3);
    }

    printf("\n******************************************************************************************\n");
    printf("* Reading slipmodel file successful:                                                     *\n");
    printf("* %-87s*\n", (*pscmp).filename);
    printf("*                                                                                        *\n");
    printf("* number of segments = %-10d number of patches  = %-10d                        *\n", (*pscmp).numberSegments, sum_no_patches);
    printf("******************************************************************************************\n");

    fclose(fin);

    strcpy((*pscmp).grnFolder, "-");


    // set patch indices
    int idum = 1;
    (*pscmp).indicesPatches = imatrix(1,(*pscmp).numberSegments,1,2);
    for (i=1; i<=(*pscmp).numberSegments; i++) {
        (*pscmp).indicesPatches[i][1] = idum;
        (*pscmp).indicesPatches[i][2] = idum + (*pscmp).numberPatches[i] - 1;
        idum = (*pscmp).indicesPatches[i][2] + 1;
    }

    (*pscmp).numberPatchesTotal = sum_no_patches;

    return (0);

}




/******************************************************************/
/*                WRITE SLIPMODEL TO FILE                         */
/******************************************************************/

int slipmodelWrite(struct slipmodel pscmp) {

    int sum_patches, patch, segment, **borders, i;
    FILE *fid;


    fid = fopen(pscmp.filename, "w");
    if (fid == NULL) {
        fprintf(stderr, "\nError opening slipmodel file %s\n(in function writeSlipmodel)\n", pscmp.filename);
        return(1);
    }

    fprintf(fid, "#===============================================================================\n");
    fprintf(fid, "# This is input file of FORTRAN77 program \"pscmp08\" for modeling post-seismic\n");
    fprintf(fid, "# deformation induced by earthquakes in multi-layered viscoelastic media using\n");
    fprintf(fid, "# the Green\'s function approach. The earthquke source is represented by an\n");
    fprintf(fid, "# arbitrary number of rectangular dislocation planes. For more details, please\n");
    fprintf(fid, "# read the accompanying READ.ME file.\n");
    fprintf(fid, "#\n");
    fprintf(fid, "# written by Rongjiang Wang\n");
    fprintf(fid, "# GeoForschungsZentrum Potsdam\n");
    fprintf(fid, "# e-mail: wang@gfz-potsdam.de\n");
    fprintf(fid, "# phone +49 331 2881209\n");
    fprintf(fid, "# fax +49 331 2881204\n");
    fprintf(fid, "#\n");
    fprintf(fid, "# Last modified: Potsdam, July, 2008\n");
    fprintf(fid, "#\n");
    fprintf(fid, "#################################################################\n");
    fprintf(fid, "##                                                             ##\n");
    fprintf(fid, "## Green\'s functions should have been prepared with the        ##\n");
    fprintf(fid, "## program \"psgrn08\" before the program \"pscmp08\" is started.  ##\n");
    fprintf(fid, "##                                                             ##\n");
    fprintf(fid, "## For local Cartesian coordinate system, the Aki\'s convention ##\n");
    fprintf(fid, "## is used, that is, x is northward, y is eastward, and z is   ##\n");
    fprintf(fid, "## downward.                                                   ##\n");
    fprintf(fid, "##                                                             ##\n");
    fprintf(fid, "## If not specified otherwise, SI Unit System is used overall! ##\n");
    fprintf(fid, "##                                                             ##\n");
    fprintf(fid, "#################################################################\n");
    fprintf(fid, "#===============================================================================\n");
    fprintf(fid, "# OBSERVATION ARRAY\n");
    fprintf(fid, "# =================\n");
    fprintf(fid, "# 1. selection for irregular observation positions (= 0) or a 1D observation\n");
    fprintf(fid, "#    profile (= 1) or a rectangular 2D observation array (= 2): iposrec\n");
    fprintf(fid, "#\n");
    fprintf(fid, "#    IF (iposrec = 0 for irregular observation positions) THEN\n");
    fprintf(fid, "#\n");
    fprintf(fid, "# 2. number of positions: nrec\n");
    fprintf(fid, "#\n");
    fprintf(fid, "# 3. coordinates of the observations: (lat(i),lon(i)), i=1,nrec\n");
    fprintf(fid, "#\n");
    fprintf(fid, "#    ELSE IF (iposrec = 1 for regular 1D observation array) THEN\n");
    fprintf(fid, "#\n");
    fprintf(fid, "# 2. number of position samples of the profile: nrec\n");
    fprintf(fid, "#\n");
    fprintf(fid, "# 3. the start and end positions: (lat1,lon1), (lat2,lon2)\n");
    fprintf(fid, "#\n");
    fprintf(fid, "#    ELSE IF (iposrec = 2 for rectanglular 2D observation array) THEN\n");
    fprintf(fid, "#\n");
    fprintf(fid, "# 2. number of x samples, start and end values: nxrec, xrec1, xrec2\n");
    fprintf(fid, "#\n");
    fprintf(fid, "# 3. number of y samples, start and end values: nyrec, yrec1, yrec2\n");
    fprintf(fid, "#\n");
    fprintf(fid, "#    sequence of the positions in output data: lat(1),lon(1); ...; lat(nx),lon(1);\n");
    fprintf(fid, "#    lat(1),lon(2); ...; lat(nx),lon(2); ...; lat(1),lon(ny); ...; lat(nx),lon(ny).\n");
    fprintf(fid, "#\n");
    fprintf(fid, "#    Note that the total number of observation positions (nrec or nxrec*nyrec)\n");
    fprintf(fid, "#    should be <= NRECMAX (see pecglob.h)!\n");
    fprintf(fid, "#===============================================================================\n");

    // OBSERVATION ARRAY
    fprintf(fid, "  2\n");
    fprintf(fid, "%-4d %.6lf  %.6lf\n", pscmp.numberObsLat, pscmp.obsLatMin, pscmp.obsLatMax);
    fprintf(fid, "%-4d %.6lf  %.6lf\n", pscmp.numberObsLon, pscmp.obsLonMin, pscmp.obsLonMax);

    fprintf(fid, "#===============================================================================\n");
    fprintf(fid, "# OUTPUTS\n");
    fprintf(fid, "# =======\n");
    fprintf(fid, "#\n");
    fprintf(fid, "# 1. select (1/0) output for los displacement (only for snapshots, see below),\n");
    fprintf(fid, "#    x, y, and z-cosines to the INSAR orbit: insar, xlos, ylos, zlos\n");
    fprintf(fid, "#\n");
    fprintf(fid, "#    if this option is selected, the snapshots will include additional data:\n");
    fprintf(fid, "#    LOS_Dsp = los displacement to the given satellite orbit.\n");
    fprintf(fid, "#\n");
    fprintf(fid, "# 2. select (1/0) output for Coulomb stress changes (only for snapshots, see\n");
    fprintf(fid, "#    below): icmb, friction, Skempton ratio, strike, dip, and rake angles [deg]\n");
    fprintf(fid, "#    describing the uniform regional master fault mechanism, the uniform regional\n");
    fprintf(fid, "#    principal stresses: sigma1, sigma2 and sigma3 [Pa] in arbitrary order (the\n");
    fprintf(fid, "#    orietation of the pre-stress field will be derived by assuming that the\n");
    fprintf(fid, "#    master fault is optimally oriented according to Coulomb failure criterion)\n");
    fprintf(fid, "#\n");
    fprintf(fid, "#    if this option is selected (icmb = 1), the snapshots will include additional\n");
    fprintf(fid, "#    data:\n");
    fprintf(fid, "#    CMB_Fix, Sig_Fix = Coulomb and normal stress changes on master fault;\n");
    fprintf(fid, "#    CMB_Op1/2, Sig_Op1/2 = Coulomb and normal stress changes on the two optimally\n");
    fprintf(fid, "#                       oriented faults;\n");
    fprintf(fid, "#    Str_Op1/2, Dip_Op1/2, Slp_Op1/2 = strike, dip and rake angles of the two\n");
    fprintf(fid, "#                       optimally oriented faults.\n");
    fprintf(fid, "#\n");
    fprintf(fid, "#    Note: the 1. optimally orieted fault is the one closest to the master fault.\n");
    fprintf(fid, "#\n");
    fprintf(fid, "# 3. output directory in char format: outdir\n");
    fprintf(fid, "#\n");
    fprintf(fid, "# 4. select outputs for displacement components (1/0 = yes/no): itout(i), i=1,3\n");
    fprintf(fid, "#\n");
    fprintf(fid, "# 5. the file names in char format for the x, y, and z components:\n");
    fprintf(fid, "#    toutfile(i), i=1,3\n");
    fprintf(fid, "#\n");
    fprintf(fid, "# 6. select outputs for stress components (1/0 = yes/no): itout(i), i=4,9\n");
    fprintf(fid, "#\n");
    fprintf(fid, "# 7. the file names in char format for the xx, yy, zz, xy, yz, and zx components:\n");
    fprintf(fid, "#    toutfile(i), i=4,9\n");
    fprintf(fid, "#\n");
    fprintf(fid, "# 8. select outputs for vertical NS and EW tilt components, block rotation, geoid\n");
    fprintf(fid, "#    and gravity changes (1/0 = yes/no): itout(i), i=10,14\n");
    fprintf(fid, "#\n");
    fprintf(fid, "# 9. the file names in char format for the NS tilt (positive if borehole top\n");
    fprintf(fid, "#    tilts to north), EW tilt (positive if borehole top tilts to east), block\n");
    fprintf(fid, "#    rotation (clockwise positive), geoid and gravity changes: toutfile(i), i=10,14\n");
    fprintf(fid, "#\n");
    fprintf(fid, "#    Note that all above outputs are time series with the time window as same\n");
    fprintf(fid, "#    as used for the Green\'s functions\n");
    fprintf(fid, "#\n");
    fprintf(fid, "#10. number of scenario outputs (\"snapshots\": spatial distribution of all above\n");
    fprintf(fid, "#    observables at given time points; <= NSCENMAX (see pscglob.h): nsc\n");
    fprintf(fid, "#\n");
    fprintf(fid, "#11. the time [day], and file name (in char format) for the 1. snapshot;\n");
    fprintf(fid, "#12. the time [day], and file name (in char format) for the 2. snapshot;\n");
    fprintf(fid, "#13. ...\n");
    fprintf(fid, "#\n");
    fprintf(fid, "#    Note that all file or directory names should not be longer than 80\n");
    fprintf(fid, "#    characters. Directories must be ended by / (unix) or \\ (dos)!\n");
    fprintf(fid, "#===============================================================================\n");

    // OUTPUTS
    fprintf(fid, " 0  0.000000  0.000000  0.000000\n");
    fprintf(fid, " 1  %-6.3lf %-6.3lf %-8.3lf %-8.3lf %-8.3lf %-8.3lG %-8.3lG %-8.3lG\n", pscmp.friction, pscmp.skempton, pscmp.strike, pscmp.dip, pscmp.rake, pscmp.sigma1, pscmp.sigma2, pscmp.sigma3);
    fprintf(fid, "\'%s\'\n", pscmp.cmbOutFolder); // output folder for cmb calculation
    fprintf(fid, "  0                 0                0\n");
    fprintf(fid, "  \'displ_x.dat\'    \'displ_y.dat\'    \'displ_z.dat\'\n");
    fprintf(fid, "  0                  0                  0                 0                  0                 0\n");
    fprintf(fid, "  \'stress_xx.dat\'   \'stress_yy.dat\'   \'stress_zz.dat\'   \'stress_xy.dat\'    \'stress_yz.dat\'   \'stress_zx.dat\'\n");
    fprintf(fid, "  0               0              0            0            0\n");
    fprintf(fid, "  \'tiltx.dat\'    \'tilty.dat\'    \'rot.dat\'   \'gd.dat\'    \'gr.dat\'\n");
    fprintf(fid, "  1\n");
    fprintf(fid, "     0.00  \'%s\'\n", pscmp.cmbOutFile); // output file name for cmb calculation
    fprintf(fid, "#===============================================================================\n");
    fprintf(fid, "#\n");
    fprintf(fid, "# GREEN\'S FUNCTION DATABASE\n");
    fprintf(fid, "# =========================\n");
    fprintf(fid, "# 1. directory where the Green\'s functions are stored: grndir\n");
    fprintf(fid, "#\n");
    fprintf(fid, "# 2. file names (without extensions!) for the 13 Green\'s functions:\n");
    fprintf(fid, "#    3 displacement komponents (uz, ur, ut): green(i), i=1,3\n");
    fprintf(fid, "#    6 stress components (szz, srr, stt, szr, srt, stz): green(i), i=4,9\n");
    fprintf(fid, "#    radial and tangential components measured by a borehole tiltmeter,\n");
    fprintf(fid, "#    rigid rotation around z-axis, geoid and gravity changes (tr, tt, rot, gd, gr):\n");
    fprintf(fid, "#    green(i), i=10,14\n");
    fprintf(fid, "#\n");
    fprintf(fid, "#    Note that all file or directory names should not be longer than 80\n");
    fprintf(fid, "#    characters. Directories must be ended by / (unix) or \\ (dos)! The\n");
    fprintf(fid, "#    extensions of the file names will be automatically considered. They\n");
    fprintf(fid, "#    are \".ep\", \".ss\", \".ds\" and \".cl\" denoting the explosion (inflation)\n");
    fprintf(fid, "#    strike-slip, the dip-slip and the compensated linear vector dipole\n");
    fprintf(fid, "#    sources, respectively.\n");
    fprintf(fid, "#\n");
    fprintf(fid, "#===============================================================================\n");
    // Green's function folder
    fprintf(fid, " \'%s\'\n", pscmp.grnFolder);
    fprintf(fid, " \'uz\'  \'ur\'  \'ut\'\n");
    fprintf(fid, " \'szz\' \'srr\' \'stt\' \'szr\' \'srt\' \'stz\'\n");
    fprintf(fid, " \'tr\'  \'tt\'  \'rot\' \'gd\'  \'gr\'\n");
    fprintf(fid, "#===============================================================================\n");
    fprintf(fid, "# RECTANGULAR SUBFAULTS\n");
    fprintf(fid, "# =====================\n");
    fprintf(fid, "# 1. number of subfaults (<= NSMAX in pscglob.h), latitude [deg] and east\n");
    fprintf(fid, "#    longitude [deg] of the regional reference point as  origin of the Cartesian\n");
    fprintf(fid, "#    coordinate system: ns, lat0, lon0\n");
    fprintf(fid, "#\n");
    fprintf(fid, "# 2. parameters for the 1. rectangular subfault: geographic coordinates\n");
    fprintf(fid, "#    (O_lat, O_lon) [deg] and O_depth [km] of the local reference point on\n");
    fprintf(fid, "#    the present fault plane, length (along strike) [km] and width (along down\n");
    fprintf(fid, "#    dip) [km], strike [deg], dip [deg], number of equi-size fault\n");
    fprintf(fid, "#    patches along the strike (np_st) and along the dip (np_di) (total number of\n");
    fprintf(fid, "#    fault patches = np_st x np_di), and the start time of the rupture; the\n");
    fprintf(fid, "#    following data lines describe the slip distribution on the present sub-\n");
    fprintf(fid, "#    fault:\n");
    fprintf(fid, "#\n");
    fprintf(fid, "#    pos_s[km]  pos_d[km]  slip_along_strike[m]  slip_along_dip[m]  opening[m]\n");
    fprintf(fid, "#\n");
    fprintf(fid, "#    where (pos_s,pos_d) defines the position of the center of each patch in\n");
    fprintf(fid, "#    the local coordinate system with the origin at the reference point:\n");
    fprintf(fid, "#    pos_s = distance along the length (positive in the strike direction)\n");
    fprintf(fid, "#    pos_d = distance along the width (positive in the down-dip direction)\n");
    fprintf(fid, "#\n");
    fprintf(fid, "#\n");
    fprintf(fid, "# 3. ... for the 2. subfault ...\n");
    fprintf(fid, "# ...\n");
    fprintf(fid, "#                   N\n");
    fprintf(fid, "#                  /\n");
    fprintf(fid, "#                 /| strike\n");
    fprintf(fid, "#                +------------------------\n");
    fprintf(fid, "#                |\\        p .            \\ W\n");
    fprintf(fid, "#                :-\\      i .              \\ i\n");
    fprintf(fid, "#                |  \\    l .                \\ d\n");
    fprintf(fid, "#                :90 \\  S .                  \\ t\n");
    fprintf(fid, "#                |-dip\\  .                    \\ h\n");
    fprintf(fid, "#                :     \\. | rake               \\\n");
    fprintf(fid, "#                Z      -------------------------\n");
    fprintf(fid, "#                              L e n g t h\n");
    fprintf(fid, "#\n");
    fprintf(fid, "#    Note that a point inflation can be simulated by three point openning\n");
    fprintf(fid, "#    faults (each causes a third part of the volume of the point inflation)\n");
    fprintf(fid, "#    with orientation orthogonal to each other. the results obtained should\n");
    fprintf(fid, "#    be multiplied by a scaling factor 3(1-nu)/(1+nu), where nu is the Poisson\n");
    fprintf(fid, "#    ratio at the source. The scaling factor is the ratio of the seismic\n");
    fprintf(fid, "#    moment (energy) of an inflation source to that of a tensile source inducing\n");
    fprintf(fid, "#    a plate openning with the same volume change.\n");
    fprintf(fid, "#===============================================================================\n");
    fprintf(fid, "# n_faults (Slip model by Ji Chen, USGS)\n");
    fprintf(fid, "#-------------------------------------------------------------------------------\n");

    // number of segments
    fprintf(fid, " %d\n", pscmp.numberSegments);

    fprintf(fid, "#-------------------------------------------------------------------------------\n");
    fprintf(fid, "# n   O_lat   O_lon    O_depth length  width strike dip   np_st np_di start_time\n");
    fprintf(fid, "# [-] [deg]   [deg]    [km]    [km]     [km] [deg]  [deg] [-]   [-]   [day]\n");
    fprintf(fid, "#     pos_s   pos_d    slp_stk slp_dip open\n");
    fprintf(fid, "#     [km]    [km]     [m]     [m]     [m]\n");
    fprintf(fid, "#-------------------------------------------------------------------------------\n");

    // write PATCH and SEGMENT data
    sum_patches = 0;
    borders = imatrix(1,pscmp.numberSegments,1,2);
    for (i = 1; i <= pscmp.numberSegments; i++) {
        sum_patches = sum_patches + pscmp.numberPatches[i];
        borders[i][1] = sum_patches - pscmp.numberPatches[i] + 1;
        borders[i][2] = sum_patches;
    }

    for (segment = 1; segment <= pscmp.numberSegments; segment++) {
        fprintf(fid, " %-3d %-12.6lf %-12.6lf %-6.2lf %-6.2lf %-6.2lf %-8.2lf %-8.2lf %-5.0lf %-5.0lf %-12.5lf\n", segment,
                pscmp.segmentData[segment][1], pscmp.segmentData[segment][2], pscmp.segmentData[segment][3], pscmp.segmentData[segment][4], pscmp.segmentData[segment][5],
                pscmp.segmentData[segment][6], pscmp.segmentData[segment][7], pscmp.segmentData[segment][8], pscmp.segmentData[segment][9], pscmp.segmentData[segment][10]);
        for (patch = borders[segment][1]; patch <= borders[segment][2]; patch++) {
            fprintf(fid, "   %10.4lf %10.4lf %10.4lf %10.4lf %10.4lf\n", pscmp.patchData[patch][1],
                    pscmp.patchData[patch][2], pscmp.patchData[patch][3], pscmp.patchData[patch][4], pscmp.patchData[patch][5]);
        }
    }

    fprintf(fid, "#================================end of input===================================\n");

    fclose(fid);

    printf("\n****************************************************************************************************\n");
    printf("* Writing pscmp file successful:                                                                   *\n");
    printf("* %-97s*\n", pscmp.filename);
    printf("*                                                                                                  *\n");
    printf("* number of segments = %-10d number of patches  = %-10d                                  *\n", pscmp.numberSegments, sum_patches);
    printf("****************************************************************************************************\n");

    free_imatrix(borders,1,pscmp.numberSegments,1,2);

    return(0);

}



/******************************************************************/
/*                     COPY SLIPMODEL                             */
/******************************************************************/

// copy the content of a slipmodel struct to another slipmodel struct

int slipmodelCopy(struct slipmodel original, struct slipmodel *copy) {

    int i, j;

    // slipmodel filename
    strcpy((*copy).filename, original.filename);

    // segment data (create and fill matrix)
    (*copy).segmentData = dmatrix(1, original.numberSegments, 1, 10);
    for (i=1; i<= original.numberSegments; i++) {
        for (j=1; j<=10; j++) {
            (*copy).segmentData[i][j] = original.segmentData[i][j];
        }
    }

    // patch data (create and fill matrix)
    (*copy).patchData   = dmatrix(1, original.indicesPatches[original.numberSegments][2], 1, 5);
    for (i=1; i<= original.indicesPatches[original.numberSegments][2]; i++) {
        for (j=1; j<=5; j++) {
            (*copy).patchData[i][j] = original.patchData[i][j];
        }
    }

    // number of segments
    (*copy).numberSegments = original.numberSegments;

    // number of patches vector
    (*copy).numberPatches = ivector(1, original.numberSegments);
    for (i=1; i<=original.numberSegments; i++) {
        (*copy).numberPatches[i] = original.numberPatches[i];
    }

    // indices of patches for every segment
    (*copy).indicesPatches = imatrix(1, original.numberSegments, 1, 2);
    for (i=1; i<=original.numberSegments; i++) {
        (*copy).indicesPatches[i][1] = original.indicesPatches[i][1];
        (*copy).indicesPatches[i][2] = original.indicesPatches[i][2];
    }

    (*copy).numberPatchesTotal = original.numberPatchesTotal;
    (*copy).numberObsLat = original.numberObsLat;
    (*copy).obsLatMin = original.obsLatMin;
    (*copy).obsLatMax = original.obsLatMax;
    (*copy).numberObsLon = original.numberObsLon;
    (*copy).obsLonMin = original.obsLonMin;
    (*copy).obsLonMax = original.obsLonMax;
    (*copy).friction = original.friction;
    (*copy).skempton = original.skempton;
    (*copy).strike = original.strike;
    (*copy).dip = original.dip;
    (*copy).rake = original.rake;
    (*copy).sigma1 = original.sigma1;
    (*copy).sigma2 = original.sigma2;
    (*copy).sigma3 = original.sigma3;

    strcpy((*copy).cmbOutFolder, original.cmbOutFolder);
    strcpy((*copy).cmbOutFile, original.cmbOutFile);
    strcpy((*copy).grnFolder, original.grnFolder);


    return(0);
}





int slipmodelFree(struct slipmodel *model) {

     free_dmatrix((*model).segmentData, 1, (*model).numberSegments, 1, 10);
     free_dmatrix((*model).patchData, 1, (*model).indicesPatches[(*model).numberSegments][2], 1, 5);
     free_ivector((*model).numberPatches, 1, (*model).numberSegments);
     free_imatrix((*model).indicesPatches, 1, (*model).numberSegments, 1, 2);

    return(0);
}



int slipmodelRandomslip(struct slipmodel *model, long seedIn, double Me, double k_exponent, int usetaper) {


    double sum_slip=0, sum_slip_orig=0, av_slip_orig, dummy, segmentAreaTotal=0, M0;

    double discx, discy, oxdim, oydim, segmentMe, *slip, *coords_x, *coords_y;

    int i, segmentNumber, no_patches_x, no_patches_y;
    long seed;

    double pi = 3.141592653589793;

    seed = seedIn;

    // estimate average slip of slip distribution
    for (i = 1; i <= (*model).numberPatchesTotal; i++) {
        sum_slip_orig += sqrt(pow((*model).patchData[i][3], 2) + pow((*model).patchData[i][4], 2));
    }
    av_slip_orig = sum_slip_orig / (*model).numberPatchesTotal;

    // total area
    for (i=1;i<=(*model).numberSegments;i++)
        segmentAreaTotal += (*model).segmentData[i][4] * (*model).segmentData[i][5];

    slip     = dvector(1, (*model).numberPatchesTotal);
    coords_x = dvector(1, (*model).numberPatchesTotal);
    coords_y = dvector(1, (*model).numberPatchesTotal);

    for (i = 1; i <= (*model).numberPatchesTotal; i++) {
        coords_x[i] = (*model).patchData[i][1];
        coords_y[i] = (*model).patchData[i][2];
    }


    for (segmentNumber = 1; segmentNumber <= (*model).numberSegments; segmentNumber++) {

        // estimate segment magnitude
        dummy = (*model).segmentData[segmentNumber][4] * (*model).segmentData[segmentNumber][5] / segmentAreaTotal;
        M0 = pow(10, 1.5 * (Me + 6.1)) * dummy;
        segmentMe = (log10(M0) / 1.5) - 6.1;
        /*******************************************************************************/

        // patch size in strike(x)/dip(y) direction
        discx = (*model).segmentData[segmentNumber][4] / (*model).segmentData[segmentNumber][8];
        discy = (*model).segmentData[segmentNumber][5] / (*model).segmentData[segmentNumber][9];

        // segment size in strike(x)/dip(y) direction
        oxdim = (*model).segmentData[segmentNumber][4];
        oydim = (*model).segmentData[segmentNumber][5];

        // number of patches in strike(x)/dip(y) direction
        no_patches_x = (*model).segmentData[segmentNumber][8];
        no_patches_y = (*model).segmentData[segmentNumber][9];



        // sort patch coordinates according to increasing x and y
        // seperated for every segment
        sort2((long) (*model).numberPatches[segmentNumber],
                &coords_x[(*model).indicesPatches[segmentNumber][1]],
                &coords_y[(*model).indicesPatches[segmentNumber][1]]);


        randslip(&slip[(*model).indicesPatches[segmentNumber][1]], &seed,
                oxdim, oydim, discx, discy,
                &coords_x[(*model).indicesPatches[segmentNumber][1]],
                &coords_y[(*model).indicesPatches[segmentNumber][1]],
                no_patches_x, no_patches_y, segmentMe, k_exponent, usetaper);

    }



    for (i = 1; i <= (*model).numberPatchesTotal; i++) {
        sum_slip += slip[i];
    }

    // normalize slip proportional to original average slip
    for (i = 1; i <= (*model).numberPatchesTotal; i++) {
        slip[i] = (slip[i] / sum_slip) * av_slip_orig * (*model).numberPatchesTotal;
    }

    // fill patchdataNew matrix
    for (i = 1; i <= (*model).numberPatchesTotal; i++) {
        (*model).patchData[i][1] = coords_x[i];
        (*model).patchData[i][2] = coords_y[i];
        (*model).patchData[i][3] = cos((pi / 180) * (*model).rake) * slip[i];
        (*model).patchData[i][4] = -sin((pi / 180) * (*model).rake) * slip[i];
    }


    // free slip vector
    free_dvector(slip, 1, (*model).numberPatchesTotal);
    free_dvector(coords_x, 1, (*model).numberPatchesTotal);
    free_dvector(coords_y, 1, (*model).numberPatchesTotal);

    return (0);

}


int randslip(double *slip, long *seed, double oxdim, double oydim, double discx, double discy, double *coords_x, double *coords_y,
        int no_patches_x, int no_patches_y, double Me, double k_exponent, int usetaper) {

    double *nx, *ny, *k_x, *k_y;

    double *REslip, *IMslip, *FTrandRE, *FTrandIM;
    double kny_x, kny_y, dk_x, dk_y, kc, kcx, kcy;

    double *specRE, *specIM, *taper;
    double minslip, alphax, alphay, sum_slip=0, dummy;

    int i, j, k, l, no_patches;
    long index;

    double pi = 3.141592653589793;

    no_patches = no_patches_x * no_patches_y;

    nx = dvector(1, no_patches);
    ny = dvector(1, no_patches);

    double min_x = dmin(no_patches, coords_x, &index);
    double min_y = dmin(no_patches, coords_y, &index);

    for (i = 1; i <= no_patches; i++) {
        nx[i] = coords_x[i-1] - min_x + 0.5 * discx;
        ny[i] = coords_y[i-1] - min_y + 0.5 * discy;
    }

    // generate wavenumber vector
    kny_x = 1 / discx / 2; // Nyquist wave numbers=maximal wave number
    kny_y = 1 / discy / 2;
    dk_x = 1 / oxdim; // discretisation in wavenumber space
    dk_y = 1 / oydim;
    k_x = dvector(1, no_patches);
    k_y = dvector(1, no_patches);
    for (i = 0; i < no_patches; i++) {
        l = floor(i / no_patches_y);
        if (l > no_patches_x / 2) {
            l = floor(i / no_patches_y) - no_patches_x;
        }
        j = i % no_patches_y;
        if (j > no_patches_y / 2) {
            j = i % no_patches_y - no_patches_y;
        }
        k_x[i + 1] = l * dk_x;
        k_y[i + 1] = j * dk_y;
    }

    //  create random field of size no_patches and make 2D-DFT of it
    REslip = dvector(1, no_patches);
    IMslip = dvector(1, no_patches);
    FTrandRE = dvector(1, no_patches);
    FTrandIM = dvector(1, no_patches);


    // initialize random generator
    dummy = ran1(&*seed);

    for (i = 1; i <= no_patches; i++) {
        REslip[i] = ran1(&*seed);
        IMslip[i] = 0;
    }
    dft2d(no_patches_x, no_patches_y, 0, REslip, IMslip, FTrandRE, FTrandIM);


    // calculate spectrum with random phases and amplitude decay 1/(k)^(2)
    // corresponding to a Hurst exponent of 1 for k higher than the corner wavenumber kc
    specRE = dvector(1, no_patches);
    specIM = dvector(1, no_patches);
    kc = pow(10, (1.82 - 0.5 * (Me))); // corner wavenumber according to Causse et al. (2010),
    kcx = 3 / oxdim;
    kcy = 3 / oydim;

    for (i = 1; i <= no_patches; i++) {
        specRE[i] = FTrandRE[i] / (pow(fabs(k_x[i]), k_exponent) + pow(fabs(k_y[i]), k_exponent));
        specIM[i] = FTrandIM[i] / (pow(fabs(k_x[i]), k_exponent) + pow(fabs(k_y[i]), k_exponent));
        if (sqrt(pow(fabs(k_x[i]), 2.0) + pow(fabs(k_y[i]), 2.0)) < kc) {
            specRE[i] = FTrandRE[i] / (pow(fabs(kc), k_exponent));
            specIM[i] = FTrandIM[i] / (pow(fabs(kc), k_exponent));
        }
        if (k_x[i] == 0 && k_y[i] == 0) {
            specRE[i] = FTrandRE[i] / (pow(fabs(kc), k_exponent));
            specIM[i] = FTrandIM[i] / (pow(fabs(kc), k_exponent));
        }
    }

    // inverse fourier transform to obtain final slip distribution
    dft2d(no_patches_x, no_patches_y, 1, specRE, specIM, REslip, IMslip);
    // set negativ values of slip to 0 by substracting the minimum value
    minslip = REslip[1];
    for (i = 1; i <= no_patches; i++) {
        if (REslip[i] < minslip) {
            minslip = REslip[i];
        }
    }

    /******************************************************************************/

    /*********/
    /* TAPER */
    /*********/
    // create taper (similar to tukey window but with onlay 1/4 of the cos period, alpha = fraction of length used for decay, alpha=1 = rectangular window, alpha=0 = similar Hann window)
    // if refdepth<1 there is no taper at the upper border
    taper = dvector(1, no_patches);
    //  printf("no_patches_x=%d, no_patches_y=%d\n", no_patches_x, no_patches_y);
    alphax = 1 / (double) no_patches_x * ((double) no_patches_x - 20);
    alphay = 1 * ((double) no_patches_y - 20) / (double) no_patches_y;
    //  printf("alphax=%lf, alphay=%lf\n", alphax, alphay);
    //if (refdepth > 1) {
    if (usetaper == 1) {
        for (i = 1; i <= no_patches; i++) {
            if (fabs(nx[i] - (no_patches_x - 1) / 2 * discx) > alphax * no_patches_x / 2 * discx && fabs(ny[i]
                    - (no_patches_y - 1) / 2 * discy) > alphay * no_patches_y / 2 * discy) {
                taper[i] = cos(pi / 2 * (fabs(fabs(nx[i] - (no_patches_x + 1) / 2 * discx) - alphax * no_patches_x / 2
                        * discx)) / ((1 - alphax) * no_patches_x / 2 * discx)) * cos(pi / 2 * (fabs(fabs(ny[i]
                        - (no_patches_y + 1) / 2 * discy) - alphay * no_patches_y / 2 * discy)) / ((1 - alphay)
                        * no_patches_y / 2 * discy));
            }
            else if (fabs(nx[i] - (no_patches_x - 1) / 2 * discx) > alphax * no_patches_x / 2 * discx && fabs(ny[i]
                    - (no_patches_y - 1) / 2 * discy) <= alphay * no_patches_y / 2 * discy) {
                taper[i] = cos(pi / 2 * (fabs(fabs(nx[i] - (no_patches_x + 1) / 2 * discx) - alphax * no_patches_x / 2
                        * discx)) / ((1 - alphax) * no_patches_x / 2 * discx));
            }
            else if (fabs(nx[i] - (no_patches_x - 1) / 2 * discx) <= alphax * no_patches_x / 2 * discx && fabs(ny[i]
                    - (no_patches_y - 1) / 2 * discy) > alphay * no_patches_y / 2 * discy) {
                taper[i] = cos(pi / 2 * (fabs(fabs(ny[i] - (no_patches_y + 1) / 2 * discy) - alphay * no_patches_y / 2
                        * discy)) / ((1 - alphay) * no_patches_y / 2 * discy));
            }
            else if (fabs(nx[i] - (no_patches_x - 1) / 2 * discx) <= alphax * no_patches_x / 2 * discx && fabs(ny[i]
                    - (no_patches_y - 1) / 2 * discy) <= alphay * no_patches_y / 2 * discy) {
                taper[i] = 1;
            }
        }
    }
    //else {
    // don't taper upper edge
    else if (usetaper == 2) {
        for (i = 1; i <= no_patches; i++) {
            if (fabs(nx[i] - (no_patches_x - 1) / 2 * discx) > alphax * no_patches_x / 2 * discx && (ny[i]
                    - (no_patches_y - 1) / 2 * discy) > alphay * no_patches_y / 2 * discy) {
                taper[i] = cos(pi / 2 * (fabs(fabs(nx[i] - (no_patches_x + 1) / 2 * discx) - alphax * no_patches_x / 2
                        * discx)) / ((1 - alphax) * no_patches_x / 2 * discx)) * cos(pi / 2 * (fabs(fabs(ny[i]
                        - (no_patches_y + 1) / 2 * discy) - alphay * no_patches_y / 2 * discy)) / ((1 - alphay)
                        * no_patches_y / 2 * discy));
            }
            else if (fabs(nx[i] - (no_patches_x - 1) / 2 * discx) > alphax * no_patches_x / 2 * discx && (ny[i]
                    - (no_patches_y - 1) / 2 * discy) <= alphay * no_patches_y / 2 * discy) {
                taper[i] = cos(pi / 2 * (fabs(fabs(nx[i] - (no_patches_x + 1) / 2 * discx) - alphax * no_patches_x / 2
                        * discx)) / ((1 - alphax) * no_patches_x / 2 * discx));
            }
            else if (fabs(nx[i] - (no_patches_x - 1) / 2 * discx) <= alphax * no_patches_x / 2 * discx && (ny[i]
                    - (no_patches_y - 1) / 2 * discy) > alphay * no_patches_y / 2 * discy) {
                taper[i] = cos(pi / 2 * (fabs(fabs(ny[i] - (no_patches_y + 1) / 2 * discy) - alphay * no_patches_y / 2
                        * discy)) / ((1 - alphay) * no_patches_y / 2 * discy));
            }
            else if (fabs(nx[i] - (no_patches_x - 1) / 2 * discx) <= alphax * no_patches_x / 2 * discx && (ny[i]
                    - (no_patches_y - 1) / 2 * discy) <= alphay * no_patches_y / 2 * discy) {
                taper[i] = 1;
            }
        }
    }

    /******************************************************************************/


    // use k-1 because pointer to first element is given which corresponds to position 0
    for (k = 1; k <= no_patches; k++) {
        if (usetaper==1 || usetaper==2)
            slip[k-1] = (REslip[k] - (minslip)) * taper[k] + 0.001; // -minslip to avoid negative slip
        else
            slip[k-1] = (REslip[k] - (minslip));
        sum_slip = sum_slip + slip[k-1];
    }

    //printf("Ohne Skalierung:\nSumme Slip = %lf   average = %lf\n", sum_slip, sum_slip / no_patches);




    free_dvector(nx, 1, no_patches);
    free_dvector(ny, 1, no_patches);

    free_dvector(k_x, 1, no_patches);
    free_dvector(k_y, 1, no_patches);

    free_dvector(REslip, 1, no_patches);
    free_dvector(IMslip, 1, no_patches);
    free_dvector(FTrandRE, 1, no_patches);
    free_dvector(FTrandIM, 1, no_patches);

    free_dvector(specRE, 1, no_patches);
    free_dvector(specIM, 1, no_patches);

    free_dvector(taper, 1, no_patches);

    return (0);

}



